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Abstract The kinematic flow pattern in slow defor¬ 
mation of a model dense granular medium is studied at 
high resolution using in situ imaging, coupled with par¬ 
ticle tracking. The deformation conflguration is inden¬ 
tation by a flat punch under macroscopic plane-strain 
conditions. Using a general analysis method, velocity 
gradients and deformation flelds are obtained from the 
disordered grain arrangement, enabling flow character¬ 
istics to be quantifled. The key observations are the 
formation of a stagnation zone, as in dilute granular 
flow past obstacles; occurrence of vortices in the flow 
immediately underneath the punch; and formation of 
distinct shear bands adjoining the stagnation zone. The 
transient and steady state stagnation zone geometry, as 
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well as the strength of the vortices and strain rates in 
the shear bands, are obtained from the experimental 
data. All of these results are well-reproduced in exact- 
scale Non-Smooth Contact Dynamics (NSCD) simula¬ 
tions. Full 3D numerical particle positions from the sim¬ 
ulations allow extraction of flow features that are ex¬ 
tremely difflcult to obtain from experiments. Three ex¬ 
amples of these, namely material free surface evolution, 
deformation of a grain column below the punch and res¬ 
olution of velocities inside the primary shear band, are 
highlighted. The variety of flow features observed in this 
model problem also illustrates the difflculty involved in 
formulating a complete micromechanical analytical de¬ 
scription of the deformation. 

PACS 83.80.Fg, 81.20.Ev, 81.05.Rm, 81.40.Lm 

Keywords dense flow, vortices, stagnation, shear 
bands, pattern formation 

1 Introduction 

Flow patterns observed in granular materials can re¬ 
semble those of fluids as well as solids, depending on the 
density [T]. At low densities, the flow is well described 
by kinetic theory [ 2113 ] . corresponding to the hydrody¬ 
namic regime. As the density is increased, long-lasting 
interparticle contacts occur, resulting in solid-like be¬ 
havior, characterized by force chains [4]. This is the 
case, for instance, in the deformation of a dense granu- 
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lar pile under the action of gravity, where the contacts 
between adjacent particles last for infinite time. 

When the particles constituting the medium are ex¬ 
tremely small, the deformation can be described by a 
suitable average over a Representative Volume Element 
(RVE). In this framework, the material is said to yield 
locally, whenever the shear stress is high [5]. However, 
as the particle size increases, the corresponding size 
of the necessary RVE also increases, and such contin¬ 
uum methods become inaccurate [6]. At the level of an 
individual particle or small clusters of particles, plas¬ 
tic flow initiation is linked to the density of collective 
structures. Examples of these are Shear Transformation 
Zones (STZs) [71[8], localized fluidization spots [9] and 
deformation/ crushing of the individual particles [5]. 

Apart from obvious interest to problems in soil me¬ 
chanics, the flow of dense granular materials is central 
to tablet processing in the pharmaceutical sector m 
and powder deformation processing of metals and ce¬ 
ramics HB. Dense granular materials have also served 
as a model for studying shear banding instabilities m 
[71l8lll3[ll4l[T5] in amorphous metals (Bulk Metallic Glasses 
or BMGs). In all of these applications, the key prob¬ 
lem is that of predicting deformation and flow patterns 
under a specified kinematic configuration. Besides en¬ 
abling improved understanding of energy dissipation 
and losses during flow, such quantitative information 
will also help in developing methods to control the level 
of homogeneity of the resulting material microstruc¬ 
ture. 

The focus of the present study is on characterizing 
dense flow in the mesoscopic regime, where the con¬ 
stituent particles are large — about 1/25^^ a typical 
macroscopic problem length — in contrast to a mate¬ 
rial such as sand, where the flow can be well approx¬ 
imated as a continuum. Eor this purpose, we use the 
model experimental system of punch indentation, with 
corresponding 1 : 1 numerical simulations. The selec¬ 
tion of this system is motivated by the varied flow fea¬ 
tures observed in punch indentation with metals as well 
as sand [5 HT6[ll7[ll8[ll9ll2Qj . The system captures many 
features common to deformation and powder process¬ 


ing. A flat-ended rectangular punch is used to penetrate 
the surface of the material, under macroscopic plane- 
strain conditions. Microscopically, however, each grain 
can move in all three dimensions. The flow is character¬ 
ized at high-resolution using in situ imaging and par¬ 
ticle tracking techniques. In order to concentrate solely 
on the role of structural rearrangements, the effects of 
particle deformation and polydispersity are suppressed 
by the use of nearly rigid monodisperse spherical grains. 
Our results indicate that the observed flow patterns are 
qualitatively unchanged even if these constraints are re¬ 
laxed. 

Goncurrent with the experiments, the use of numer¬ 
ical simulations allows us to probe kinematic details 
that are not easily accessed in the imaging setup. To 
this end, we use the Non Smooth Contact Dynamics 
(NSCD) method [^1221123] to handle the long lasting 
interparticle contacts that occur in dense flow. Mim¬ 
icking the experiments, simulations are performed with 
monodisperse rigid spherical particles. The laws gov¬ 
erning individual particle dynamics are simple and the 
observed flow patterns result solely from collective mo¬ 
tion of the grains. Due to fundamental differences be¬ 
tween granular packings in two and three dimensions 
[24] . 3D simulations are necessary to handle monodis- 
persity without causing crystallization. 

With this as background, we describe the experi¬ 
mental setup and numerical simulation method in Sec.[^ 
Our main results are presented in Sec. The primary 
flow features are then described and analyzed; these 
are quantitatively reproduced by the numerical simula¬ 
tions. In Sec.|^ we discuss the origin and characteristics 
of the flow patterns and draw comparisons with similar 
features seen in other experiments, notably in punch in¬ 
dentation of metals. The results and their implications 
are summarized in Sec. [H 

2 Methods 

In order to quantify effects of structural rearrangements 
on the kinematic flow pattern in a dense non-cohesive 
granular medium, model experiments were carried out 
using plane-strain punch indentation. Gorresponding 
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numerical simulations were performed using the Non- 
Smooth Contact Dynamics (NSCD) method. 


2.1 Experimental setup 

Figure shows a schematic of the plane-strain inden¬ 
tation, wherein a flat-ended rectangular punch, made 
of stainless steel with dimensions 35 mm x 53 mm x 
25 mm {2D x Lp x Tp), slowly penetrates the granular 
medium. The granular material consists of monodis- 
perse spherical steel balls with diameter d = 1 mm. 
These were chosen for ease of tracking and to avoid the 
influence of size variation and particle deformation on 
the observed flow pattern. The balls were placed in a 
rectangular A1 alloy box of dimensions 140 mm x 65 
mm X 25 mm, as shown in the figure. The front face was 
covered with a glass plate to enable direct observation 
and imaging. 

The thickness Tp of the punch (in the ^-direction) 
matched that of the material, confining the punch mo¬ 
tion to the xy plane. The tolerance ATp in the thick¬ 
ness dimension was maintained at ±0.1 mm {ATp d), 
thereby ensuring that particles did not get into the gap 
between the punch wall and the front face of the box. 
The origin of coordinates was fixed at the midpoint of 
the bottom (fiat) face of the punch, with the x-axis co¬ 
inciding with the bottom face, i.e. along the width 2D 
of the punch — see Fig. 0 (right). 

The deforming granular medium was imaged using a 
high-speed CMOS camera (PCO dimax), with a spatial 
resolution of 60 ym per pixel. The field of view covered 
an area of 120 mm x 50 mm, and was illuminated by a 
120 W halogen lamp light source. The experiments were 
done on a Universal Testing Machine (UTM) under ve¬ 
locity control. Three indentation velocities Vp were used 
— 0.02 mm/s, 0.2 mm/s and 2.0 mm/s — covering 3 
different orders of magnitude. Quantitative information 
about the flow field was obtained from the recorded 
images, using image analysis techniques. This involved 
first identifying the particles composing the medium 
in each frame, and then tracking them using Particle 
Tracking Velocimetry (PTV) methods. Post processing 


of experimental data is described in Sec. 2.3 


Constant velocity Vp 



(a) 



(b) 

Fig. 1 Problem setup, (a) (Left) Schematic of the experi¬ 
mental setup with sample image frame superimposed over 
the field of view (highlighted in red). (Right) Reference frame 
with punch stationary, showing the coordinate system used. 
The punch width 2D, used for normalizing lengths, is also de¬ 
picted. (b) Perspective view of the container with the granular 
material. OO' is the plane of observation. Numerical simula¬ 
tions are performed for the entire 3D domain, including the 
front and back walls, corresponding to physical dimensions. 
Simulation results are visualized only for the center plane SS'. 


2.2 Numerical simulation 

Discrete Element Method (DEM) simulation of the punch 
indentation was performed to calculate the 3D kine¬ 
matic properties of the medium. Several variants of 
DEM exist, each trading conceptual simplicity for sta¬ 
bility and timestep length [2311^1^ . The dense gran¬ 
ular material used in the experiments consisted of hun¬ 
dreds of thousands of persistent contacts. For this rea¬ 
son, we used the Non-Smooth Contact Dynamics (NSCD) 
method [mna with the open source solfec code [27] to 
carry out the simulations. 

The generalized coordinates (positions and orienta¬ 
tions) and velocities (linear and angular) of all particles 
are assembled in the vectors q, u. The complete equa- 
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tions for time-iteration can be represented as follows 
h . 


a. f 
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( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 


with the superscripts denoting the values of the vari¬ 
ables at timestep t, t + | and t + h, with step size h. 
Here M is an inertia operator, f is a vector containing 
all external forces (such as gravity) and A is the vector 
of reaction forces due to contact between particles. H 
is an operator that translates between the generalized 
coordinates and local contact coordinates [HEi]. 

The basic idea behind the NSCD scheme is two¬ 
fold. Firstly, timestepping is done from t to t + h/2, as 
in Eq. |2.1| Secondly, the coupled rigid-body constraint 
equations of the form 


C(Hu*+'*,yl) = 0 


(2.4) 


are solved for A. The constraint equation arises from 
both impenetrability of particles and static/ dynamic 
friction (Signorini and Coulomb constraints respectively). 
Eq. |2.4| is solved after constraints are updated and new 
contact points detected, using the new positions q^+^/2^ 
following Eq. |2.1[ Timestepping is then continued, by 
updating the velocities in Eq. |2.2| and positions from 


t + h/2tot + hin Eq. 2.3 This scheme enables the res¬ 
olution of a large rigid-body system via a stable and ef¬ 
ficient implicit solution, thus allowing large timesteps h. 
This approach is what differentiates the NSCD method 
from other DEM variants : it does not assume a Hertzian 
or spring contact model [23] to compute reaction forces. 
Eor more details on the implementation, see Ref. m- 
Dimensions in the simulations were chosen to cor¬ 
respond 1 : 1 with the experimental setup, shown in 
Eig. |l(b)| The simulated granular medium was com¬ 
posed of 233,069 rigid spheres. The initial packing was 
formed using a collision driven packing algorithm [29] . 
which generated a disordered packing of monodisperse 
spheres inside a cube of unit dimensions. Three such 
packings were assembled into a 3D box, with dimen¬ 
sions as in the experiments. The faces of the box were 


composed of rough rigid planes. This final ensemble was 
then allowed to relax under gravity using the NSCD 
method. A punch, also composed of rough rigid planar 
faces was then placed at the center of the relaxed sam¬ 
ple and lowered at a constant velocity, maintained after 
every time step. This was done to mimic the velocity 
control used in the experiments. 

The 1 : 1 simulation-to-experiment correspondence 
prevented a naiVe translation of physical problem di¬ 
mensions, which required very small timesteps, thereby 
increasing computation time. Hence, length and time 
were rescaled in the simulations, to ensure numerical 
stability during timestepping. Physical length Iq and 
time to were rescaled as I = and t = yto respec¬ 
tively with ^ = 10^, X = 31.645. This ensured that the 
acceleration due to gravity was g = 10//t^ = 10/o/to- 
The simulation punch velocity v = Cuo with C = 31.62. 
The coefficient of friction g values for all contacting 
surfaces — ball-ball, ball-punch and ball-wall — were 
set at /i = 0.4, typical for steel contacts. The visualized 
results are presented for the midplane of the simulation 
domain, see Eig. |l(b)[ 

2.3 Post-processing 

The experimental data consisted of a sequence of high- 
resolution images of the medium flowing past the punch. 
The punch velocity Vp was kept constant for each ex¬ 
periment, and the images were corrected to keep the 
punch stationary — see Eig. |l(a)| (right). The grayscale 
images from the camera were represented by 2D float¬ 
ing point arrays, containing intensity data in the range 
[0,1]. Each frame was then processed to obtain the po¬ 
sitions of the sphere centers and their displacements be¬ 
tween successive frames, as described in Ref. m- Since 
the frame rate was kept constant, velocities were ob¬ 
tained from inter-frame particle displacements over 15 
frames. Eor particle i, the x and y components of the 
velocity were Ui and i = 1, 2, • • • for N tracked 
particles. The positions of the particles themselves were 
(xi, yi). The scalars Ui and Vi represented two fields on a 
2D plane, defined on the set of points with coordinates 
ixi,yi). 
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Fig. 2 Calculation of incremental strain measures from posi¬ 
tion and displacement data. The particle center positions and 
inter-frame displacements are stored (I). The position infor¬ 
mation {xi,yi) is used to create a Delaunay triangulation of 
the entire ensemble (II). The 2D triangulation is converted 
to a 3D mesh by concatenating the particle’s displacement 
(shown here for u, a similar procedure is used for v) to its co¬ 
ordinates (III). The resulting 3D mesh is used to approximate 
the X and y gradients of u at each particle location. 


In the simulations, the positions {xi^yi^Zi) and ve¬ 
locities (ui^Vi^Wi) of all the particles were recorded af¬ 
ter every timestep. Deformation measures were then 
computed as necessary. Renders of particle positions 
were done using the Open Visualization Tool m- 


2. The triangulation, with the scalar field Ui deter¬ 
mined a 3D heightfield mesh with vertices at (x^, yi^Ui 
corresponding to each tracked particle i. Implicitly, 
in 3D, this heightfield surface was given hy F{x^y) = 

^ - u{x,y) = 0 (Fig. |2] (III)). 

3. The vertex normals were calculated for each trian¬ 
gle of the 3D heightfield mesh, using methods of 
standard computational geometry [32]. They were 
normalized to have unit norm and were of the form 
(nf,nf,nf). 

4. The heightfield, being of the form F{x,y) = 2 ) — 
u{x^y) = 0, had normals VF{x^y^z)^ given by 

VF(x, y, z) = {—dujdx, —dufdy, 1) (2.5) 

So for discrete points i, {du/dx)i = —nf/n| and 
{du/dy)i = —rv(/nl (when nf ^ 0) approximated 
Vtx at the tracked particle i. 


The velocity gradients Vu and Vt; at each point i were 
then used to approximate continuum rate measures for 
the medium. 


The vorticity uji at each particle location, was cal¬ 
culated as: 


uJi = 


/ du\ / dv\ 


( 2 . 6 ) 


2.3.1 Deformation rate measures 

In order to determine various kinematic flow proper¬ 
ties, an evaluation of the 2D velocity gradients Vu = 
{duldx^dujdy) and Vv = {dvjdx^dvjdy) was neces¬ 
sary. For a true continuum, these are usually approxi¬ 
mated using finite differences on a uniform grid. How¬ 
ever, given the discrete particulate nature of the medium, 
a true gradient (in the continuum sense) is not well- 
defined for a disordered grain network. Instead, velocity 
gradients were approximated as follows (See Fig.[^ the 
discussion is for Vtx; the same applies to Vv). 

1. A 2D Delaunay traingulation was performed with 
the obtained particle coordinates (xi^yi) (Fig. 
parts (I) and (II)) 


The strain rate tensor components were approxi¬ 
mated by 



(2.7) 


From these relations, the 2D strain rate invariant 
was calculated for each particle i for a particular frame 
as 




2 

3 



- Cyy) + ( 



+ ^^ly 


( 2 . 8 ) 


This provides a measure of the incremental shear 
strain in the medium. 
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2.3.2 Accumulated strain 


The strain rate tensor and vorticity were calculated 
from gradients of the instantaneous velocities. However, 
finite particle displacements during deformation inval¬ 
idated the small displacement gradient approximation 
for the strain tensor. Due to temporal averaging used to 
calculate the particle displacements, integration of the 
incremental strain rate blurred the details of the strain 
field. Instead, the accumulated strain in the material 
was computed as follows (for more details, see Refs. [H 

El). 

For each particle i, the distance IT from its near¬ 
est neighbors j was calculated after every timestep t. 
The deformation gradient for particle i was then de¬ 
termined by minimizing ~ being the 

separation in the initial configuration. The Lagrangian 
strain tensor was defined as = (F^F^ — I)/2. While 
the particles themselves were rigid, measured mate¬ 
rial deformation resulting from relative motion of neigh¬ 
boring particles. 

From E, the second shear strain invariant rji (in 3D) 
was calculated as 


Vi = 


^yz + + ^ly~ 


(2.9) 


{Eyy ~ EzzY + {^xx ~ ^zzY + {^yy ~ ^xxY^ ^ 


( 2 . 10 ) 


and provided a scalar measure of the material shear 
strain, corresponding to the particle i. It has been shown 
that r]i captures strain due to local structural rearrange¬ 
ments HUH]. For the 2D case, E^z = Eyz = Ezz = 0 in 
the relation for rn above. 

The expressions for ef and rji in terms of compo¬ 
nents of the strain rate (ea6)and strain tensors (Eab) 
are identical but their methods of computation are dif¬ 
ferent. 


3 Results 

The individual particle trajectories {xi{t)^yi{t)), veloc¬ 
ities (ui^Vi) and gradients (Vi^)^, (Vv)i (in the plane xy 



(a) 



(b) 


Fig. 3 Validity of plane-strain assumption at the mesoscopic 
scale, (a) Plot of particle velocity v-^ normal to observa¬ 
tion plane in a simulation run. This velocity decreases quite 
rapidly with distance from the front and back faces of the 
container box. (b) Colored particle paths from successive ex¬ 
perimental images. Punch velocity Vp = 0.2 mm/s. 


of observation) were extracted from the high-resolution 
image sequence (c/. Sec. 2.3). Simultaneously, 3D par¬ 
ticle positions and velocities were obtained from 1 : 1 
numerical NS CD simulations. Movie Ml (see supple¬ 
mental material [33] ) shows a typical sequence of images 
from an experiment and the corresponding simulation. 


3.1 Particle tracking in plane-strain 

Macroscopically, the punch is constrained to move in 
the x^-plane. However, microscopically, the particles 
composing the granular material are free to move in 
all three directions x, y and For media such as sand, 
where the particle size is much smaller than a typical 
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macroscopic length scale, the plane strain assumption 
suffices to describe the flow field accurately [5] . An anal¬ 
ysis of the particle velocities obtained from the simula¬ 
tions was used to verify the applicability of the plane 
flow approximation in the current experiment. 

For this purpose, particle velocities \vz\^ see 

coordinate system in Fig.[^, perpendicular to the view¬ 
ing {xy) plane, are extracted from the NSCD simula¬ 
tions, for parallel planes with normal along the z-axis. 


Figure 3(a) shows the variation of the average v-^jvp 
in each such plane, for increasing punch penetration 
depth Y. The graph is quantitatively similar for the 
different Vp considered. The average value of jvp is 
very small at the two peaks 3.5%), which occur just 
inside the front and back faces of the sample. More¬ 
over, after a short initial rise (between Y = 1.58 mm 
and Y = 3.16 mm), remains unchanged. The simu¬ 
lations thus indicate that the net contribution to flow 
in the 2 ;-direction, perpendicular to the viewing plane, 
is very small. This justifies the use of particle tracking 
methods restricted to the xy plane. 

A typical set of tracks, obtained from an analysis 
of the experimental data for Vp = 0.2 mm/s, is dis¬ 
played in Fig. 3(b)[ Each color line represents the tra¬ 
jectory of the center of an individual particle in the 
medium, over 30 successive frames. Particle tracking 
provided (xi^yi) and (ui^Vi), corresponding to particle 
center coordinates and inter-frame displacements (i.e. 
instantaneous velocities) respectively. The velocity gra¬ 
dients were approximated from this data (Sec. |2.3[). 


3.2 Stagnation zone and its evolution 

Keeping the punch stationary, 100 experimental images 
are stacked on top of each other to reveal the flow of 
individual grains, as in Fig. |4(a)| (top). Likewise, equiva¬ 
lent NSCD simulation data is rendered into images cor¬ 
responding to successive timesteps, and the images then 


stacked to obtain the bottom image shown in Fig. 4(a) 


The two images are very similar, consistent with the 
physical conditions being the same. This can also be 
seen from Movie Ml [33], which provides a visual com¬ 
parison between experiment and simulation. The grains 
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Fig. 4 Stagnation (or dead material) zone and its evolution, 
(a) Sequence of 100 experimental (top) and simulation (bot¬ 
tom) frames superimposed, showing the flow pattern. The 
streaks correspond to particle pathlines. The stagnation zone 
(white dashed line) appears motionless, as revealed by lack 
of streaklines, (b) Variation of base angle (6) of the stagna¬ 
tion zone with punch penetration depth Y. The dashed line 
represents tan“^(/i) used in the simulations. 


flow around the punch, creating tracks — particle path- 
lines — clearly seen in both composite images. 

A characteristic feature of the flow is the occurrence 
of a stagnation zone, or region of ‘dead material’, un¬ 
derneath the punch. This region consists of particles 
that are stationary relative to the punch. The zone is 
roughly triangular in shape, as seen in Fig.[^ and sym¬ 
metrical with respect to the ^-axis, consistent with the 
loading condition. The angle 0 of the base of the tri¬ 
angle provides a measure of the shape of the zone, see 
Fig.g a)| In order to measure the evolution of differ¬ 
ent composite images (indexed by k) are generated by 
combining 100 images around a particular value of Y. 
The top left and top right corners as well as the low- 
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est point of the stagnation zone are determined for each 
composite image. These correspond to the centers of the 
leftmost, rightmost and bottom-most stationary parti¬ 
cles in the composite image respectively. Two lines are 
then drawn joining the lowest point with the top cor¬ 
ners; the angle subtended is Ok^ Ok values are averaged 
(/c = 10) to give 0 for the particular Y value in Fig. 

0 changes with the punch penetration depth T, but 
the stagnation zone remains nearly triangular and sym¬ 
metric. A plot of 0 vs. T, as measured from experi¬ 
mental data, is shown in Fig. |4(b)| Also shown is the 
corresponding 0 value obtained from numerical simula¬ 
tions, for Vp = 0.5 mm/s. The simulations accurately 
reproduce the experimental stagnation zone geometry. 
The initial value 0 ^ 50° was found to be sensitive 
to the packing details, as confirmed by observing data 
from repeated simulation runs. After indentation to a 
sufficient depth, 0 tends to a limiting value tan“^ /i, 
independent of the sliding velocity. This behavior is also 
reproduced by the NSCD simulations. Here /i = 0.4 is 
the friction coefficient between all the surfaces in con¬ 
tact. We note that there was no observable flow along 
the bottom wall of the container for the Y values re¬ 
ported in Fig. |4(b)[ 

On a longer timescale, the size of the stagnation 
zone decreases slowly with T, even though 0 reaches a 
steady value. This is because friction causes a layer of 
particles immediately adjoining the face of the punch 
to move sideways. However, on the scale of a few (^5) 
mm, the change in the stagnation zone width is negli¬ 
gible compared to T, hence the particles do not form 
distinct pathlines in the composite image. 

3.3 Rotation field and vorticity 

In the coordinate system fixed to the punch, flow sepa¬ 
ration occurs at the bottom face of the punch. The vor¬ 
ticity uji (Eq. |2.6[ ) was used to quantify the net rotation 
of the medium, induced by differential displacements of 
adjacent particles. 

Figure shows a plot of the vorticity, linearly inter¬ 
polated from uji values at each particle location i. uji is 
scaled to lie between -1 (blue, clockwise) and +1 (red. 


Expt. 



-1.0 (CW) 


1.0 (CON) 


(a) 



(b) 


Fig. 5 Particle velocities and vorticity field near the punch. 
Arrows represent velocity of individual particles, estimated 
from successive images; arrow size is proportional to veloc¬ 
ity magn itude. Background is colored based on vorticity uoi 
(Eq. |2.6| ) using linear interpolation from values at particle 
locations, uoi is scaled between -1 (blue, clockwise) and 1 
(red, counterclockwise), (a) Experimental image for Vp = 0.2 
mm/s. Arrows point to the vortex centers, while stagnation 
zone is demarcated by black dashed lines. The layers used 
for vorticity plot are also indicated, (b) Corresponding image 
from NSCD simulation — both the stagnation zone geome¬ 
try and vortex centers are quantitatively reproduced in the 
simulations. 



x/D 

Fig. 6 Scatter plot of rotation field approximated from ex¬ 
periments, in three horizontal layers adjoining the bottom 
face of the punch. The two peaks in the layer just below the 
punch show the occurence of vortices at x/D ~ ±1.25. 
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counterclockwise). Individual particle velocities are su¬ 
perimposed on top using arrows. The arrow size reflects 
the velocity magnitude. Figure . [5(a)[ shows the result for 
an experiment with Vp = 0.2 mm/s. Qualitatively simi¬ 
lar plots were obtained for the other Vp. For a compar¬ 
ison, an equivalent plot using simulation data is shown 
in Fig. |5(b)| In this latter case, particles in a narrow 
region (of width 1.5 particle diameters) were isolated 
at the center along the 2 ;-direction (c/. Figs. |l(b) and 
1 ^ - 


The stagnation zone is marked by dashed lines in 
Fig.|5( a)| Far away from the punch, the material flows 
with velocity Up, as indicated in the figure. Two cen¬ 
ters of vorticity are observed, indicated by the two ar¬ 
rows, corresponding to regions with high uji. These two 
vortices are of opposite orientation, due to flow sepa¬ 
ration at the punch bottom face. As the material ap¬ 
proaches the bottom face, individual particles are forced 
to move away from the centerline (x = 0), thereby flow¬ 
ing around the side faces of the punch. 


To better understand the variation in material ro¬ 
tation along the y direction, uji from the experimen¬ 
tal data in three neighboring spatially fixed layers (see 
Fig. |5(a)| ) is plotted in Fig.[^ These layers are of equal 
width and adjacent to the bottom face of the punch. 
The scatter values of uJi in a layer immediately above 
the bottom face (Layer 1) are plotted in green. The val¬ 
ues in a layer right below the punch (Layer 2) are shown 
in red. uJi in the adjacent layer (Layer 3) are plotted in 
blue. In Layer 1, a sharp increase is seen in just 
around the lateral edges of the punch. This occurs at 
x/D ±1.2. In Layer 2, two peaks, corresponding to 
the vortices in Fig. |5(a) , are observed on either side of 
the punch. These peaks occur at x/D ±1.25, imme¬ 
diately below the punch bottom face. In Layer 3, the 
vorticity is significantly diminished. Hence, finite vor¬ 
ticity remains confined to a thin layer near the bottom 
face of the punch, decaying rapidly away from it. 


The observed rotation field was also well reproduced 
in the NSCD simulations (Fig. |5(b')] ), including the two 
vortices on either side of the punch. As a result of in¬ 
terparticle friction, velocities differ substantially from 


Vp only in a region near the punch. This is the cause 
for the rapid decrease in vorticity along the ^-axis, in 
Fig. Since the simulation data were extracted in the 
middle of the container, the vortex lines are parallel to 
the 2 ;-axis and extend throughout the thickness of the 
sample. Furthermore, it was observed that the vortices 
(in both experiments and simulations) were sustained 
even after the angle 0 of the dead material zone reached 
its steady value (c/. Fig. 4(b)[ ). The vortices hence exist 
during steady flow — similar observations of vortices 
have been made in experiments with dry sand [34] . 


3.4 Shear band formation 


It is well known that fine-grained granular media un¬ 
dergo shear as they flow around a punch [51135]. Due 
to local rearrangements of particles, shear strain tends 
to get concentrated in the form of bands that demar¬ 
cate regions with a sharp difference in the velocity of 
adjacent particles. 

An estimate of velocity gradients resulting in shear 


was obtained from the invariant measure ef (Eq. 2.8). 


The ef field, linearly interpolated between particle lo¬ 
cations, along with superimposed particle velocities, is 
shown in Fig. |7(a)[ The two white arrows indicate re¬ 
gions of high shear rate, occuring right below the punch. 
These regions are seen to coincide with a sharp gradient 
in particle velocity magnitude. Moreover, the fluctua¬ 
tions in ef, evident in Fig. 7(a), reflect local particle 
rearrangements. 

A scatter plot of ef, corresponding to the three lay¬ 
ers depicted in Fig. |7(a)] is shown in Fig. |7(b)] In Layer 
1, just above the bottom face of the punch (green), large 
shear strains are observed very close to the sides of the 
punch. In Layer 2 (red), high values of ef occur below 
the punch, due to particle rearrangements. Finally, sig¬ 
nificant shear rates are also seen in Layer 3 (blue). This 
is in contrast to the vorticity, which decreases rapidly 
in this layer (c/. Fig. |^. 

Regions of high velocity gradient occur at nearly 
constant spatial locations, causing shear strain accu¬ 
mulation in the form of continuously evolving narrow 
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Fig. 7 Material strain rate field from experimental data, (a) 
Half the (symmetric) flow field near the punch with particle 
velocities (a rrow s) superimposed on shear strain rate invari¬ 
ant ef (Eq. 2.8). (b) Scatter plot of ef in the three layers 
shown in (a). High strain rates accompany large changes in 
neighboring particle velocities. 


bands. A quantitative measure of this strain is pro¬ 
vided by the strain invariant r]i (Eq. |2.9[ ) . A rendered 
sequence from an NSCD simulation, with particles col¬ 
ored by their rji value, is shown in Fig. Movie M2 
[33] shows a time sequence of the evolution of r]i in the 
material. Each successive frame in Fig. corresponds 
to increasing punch penetration Y. Shear in the mate- 




Fig. 8 Calculated strain data from NSCD simulation for 
Vp = 0.5 mm/s; the particles are colored by the value of r]i. 
Only half of the flow field is shown, (a) Y = 2.37 mm. Strain 
accumulation begins at the two bottom edges Ai and A 2 of 
the punch, (b) Y = 5.53 mm. A shear band B is seen form¬ 
ing underneath the stagnation zone, (c) Y = 8.7 mm. The 
stagnation zone PQR is clearly seen, surrounded by regions 
of shear in two bands. Additional shear bands Ci and C 2 are 
beginning to develop below the stagnation zone. 


they run upto the bottom face of the punch. With in¬ 


rial first occurs at the two edges Ai and A 2 , near the 
bottom of the punch (Fig. |8(a)'] ). As the material contin¬ 
ues to flow, the formation of two distinct shear bands 
{B) is seen in Fig. |8(b)| The two bands are symmet¬ 
ric with respect to the vertical, and bound the incipi¬ 
ent stagnation zone underneath the punch. The ends of 
these bands coincide with the edges Ai and A 2 , so that 


creasing T, shear strain in these bands further increases 
(Fig. |8(c)'] ). This causes additional shear bands to form 
underneath the stagnation zone, see points Ci and C 2 
in Fig. 8(c)| At this stage, the stagnation zone geome¬ 
try is near its steady state value (i.e. 0 nearly constant). 
The total strain in this zone is much lesser than in the 
material immediately surrounding it, even though small 
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(a) 



(b) 

Fig. 9 Build-up of shear in a narrow region, (a) Spatial lo¬ 
cation for tracking particle motion. The plane AA'BB' passes 
through the thickness of the sample (in the ^-direction), (b) 
Material in the plane AA'BB', viewed along the plane normal. 
Arrows represent in-plane particle velocity i’ll. The color de¬ 
picts linearly interpolated strain invariant rj for the medium. 
Direction of particle motion indicates reduction in density 
along the band. 


local fluctuations in strain rate are observed — as seen 
in Movie M2 and Fig. 

These observations suggest that large shear strains 
are concomitant with steep velocity gradients in the 
flow. Further evidence is provided by observing par¬ 
ticles in the plane of the shear band itself. To this 
end, particles at the spatial location of the band B 
(Fig. |8(b)| , marked by the plane AA'BB' in Fig. 9(a)[ 
are isolated. This plane is chosen to subtend an angle 
(j) = arctan(4/3) with the x-axis. When viewed along 
the normal to this plane, the constituent particles form 
a thin rectangular layer, with width AA' equal to the 


thickness of the granular material. The in-plane veloc¬ 
ity of each particle i is calculated from the x and y 
velocities v^^Vy as 

= Vxcoscj)Vysmcf) (3.1) 


The velocities in the plane AA'BB' are shown by ar¬ 
rows in Fig. 9(b)[ with Y increasing from left to right. 
The background color in each of these images is calcu¬ 
lated from the strain invariant rj for the medium, ob¬ 
tained using linear interpolation from the particles’ in¬ 
dividual rji values. Initially, for small T, a zone of shear 
develops near the bottom edges of the punch, close to 
the end AA'. This corresponds to the material near A 2 


in Fig. 8(a) As the plane AA'BB' was chosen to be 
fixed in space, this zone of shear moves downward with 
increasing Y. The cause for this is the material pileup 
on either side of the punch (see sequence in Fig.|^. The 
zone of shear further increases in size with increasing 
T, seen by its extent for Y = 3.16 mm and Y = 7.90 
mm in Fig. |9(b)| 

Another interesting feature observed in Fig. 9(b)| is 
the particle velocity distribution. While particles below 
the shear zone remain stationary, those above the zone 
continue moving upward; cf. arrow sizes in Fig. |9(b)| 
for T = 6.32 mm and Y = 7.90 mm. The shear zone 
hence corresponds to high velocity gradients in the ma¬ 
terial, and low shear exists even when the particles in 
the material move with significant velocity. The mo¬ 
tion of particles away from the end BB' of the plane 
AA'BB' also results in significant local dilatancy in the 
material. Consequently, highly sheared regions result 
from irreversible motion of the consituent particles, in 
turn causing volume change in the material. 


3.5 Other structural changes 

In addition to the flow features discussed above, dense 
flow in the punch indentation configuration shows other 
interesting structural changes. Since these changes can¬ 
not be observed directly by 2D imaging, they are ex¬ 
tracted from the 3D simulation data. Two such features 
are highlighted and discussed. 
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x/D 


Fig. 10 Initial {Y = 0) and final {Y = 14.32 mm) free surface 
profile of the medium. New surface formation along the sides 
of the punch occurs by movement of grains from the bulk. 
Dotted line denotes final free surface position in the absence 
of the punch. 


3.5.1 Free surfaee evolution 

In punch indentation of metals, it is known that new 
surfaces are created along the sides of the punch, with 
increasing penetration Y [36]. This corresponds to the 
‘cutting’ mode of deformation, as opposed to the ‘ra¬ 
dial compression’ mode of deformation, more typical 
with blunt wedges. This free surface formation is ac¬ 
companied by a pileup of material along the sides of 
the punch. 

Likewise, during deformation and flow of the gran¬ 
ular medium, the free surface evolves, with a similar 
pileup around the sides of the punch (see the sequence 
in Fig.[^. In order to study the evolution of the free sur¬ 
face, particles in a thin layer (of width l.bd) adjoining 
the x-axis {y = 0) are identifled at t = 0 and followed 
during flow. The initial and flnal positions (punch pen¬ 
etration T = 0 mm and Y = 14.32 mm respectively) 
of these particles are shown as scatter plots in Fig. 

The particles initially form a horizontal layer, shown as 
black/ green triangles, in the figure. The flnal particle 
locations, shown by red/ blue circles in the figure, fol¬ 
low the shape of the punch, while forming an upward 
arch beyond x/D = ±1. The variation in thickness of 
this flnal layer reveals gaps in the surface adjoining ei¬ 
ther side of the punch. These gaps do not exist in the 





Y = 0nnm Y = 7.11 mm Y= 15.02 mm 

Fig. 11 Change in shape of vertical column of particles un¬ 
derneath the punch. The surrounding material is not shown. 
The column appears to ‘buckle’ before flowing past the sides 
of the punch. 

final material, as they are filled by particles from the 
bulk. Therefore, similar to metals, an equivalent cutting 
mode of deformation also occurs in granular materials. 

The height of the material pileup on the sides of 
the punch is measured from the hypothetical free sur¬ 
face position, shown by a dotted line in Fig. This 
is where the free surface would be if the material did 
not encounter the punch. At x/D = ±4, the two levels 
nearly coincide, confirming negligible flow against the 
side walls. 


3.5.2 Column buekling and flow 


The material immediately underneath the punch is com¬ 
pressed, being driven against its bottom face. To study 
this confined compression, a vertical column of particles 
under the punch is isolated. These particles are identi¬ 
fied at T = 0 and their subsequent positions recorded. 


a rendering of which is shown in Fig. 11 The particles 
surrounding this column are not shown in the figure. 
Due to compression by the punch, the initially vertical 
(at T = 0 mm) column develops a bulge (at T = 7.11 
mm in the figure) along its sides. With further penetra¬ 
tion, the extent of the bulge increases, until material at 
the edges finally begins to flow along the sides of the 
punch (T = 15.02 mm). Beyond this stage, material 
particles along the punch sides continue to flow upward. 
These particles hence filled the gaps that appeared in 


the ‘top surface’ of the material (see Sec. 3.5.1). Thus, 
particles that constitute the piled up material originate 
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both from the bulk away from the punch as well as from 
immediately below the punch. 

The width of the initial {Y = 0) vertical column 
in Fig. was chosen to coincide with the final width 
of the stagnation zone underneath the punch. However, 
it was observed that this choice was not crucial to the 
behavior shown in the figure — changing the initial 
column width only altered the Y at which flow past the 
sides of the punch occurred. 

4 Discussion 

Our experiments on a model granular medium have 
shown a number of unique flow features during defor¬ 
mation and plastic flow. Quantitative details of these 
features were extracted from high resolution in situ ob¬ 
servations, coupled with numerical simulations. Char¬ 
acteristics that were not amenable to analysis by 2D 
imaging were extracted from the 3D NS CD simulations. 

When viewed in the rest frame of the punch, plane- 
strain indentation resembles flow of a dense granular 
material past a static obstacle. However, an important 
difference pertains to the role of gravity — conventional 
studies of gravity-driven flow past obstacles [371[3il39] 
have involved free fall of the constituent particles. This 
reduces the average life of inter-particle contacts. In 
contrast, the flow in punch indentation occurs with mul¬ 
tiple persistent contacts. This is because the material 
flow direction is against gravity. 

The primary physical constraints in our work con¬ 
cerned the use of monodisperse, rigid spheres, in order 
to eliminate the effects of particle size variation and 
deformation. This motivated the choice for a granu¬ 
lar medium comprised of steel balls. The reported fea¬ 
tures were found to be qualitatively unchanged upon 
removing these two constraints, as seen from Movie 
M3 (supplemental material [33]), which shows punch 
indentation of a medium composed of mustard seeds. 
These are crushable polydisperse grains, yet the ob¬ 
served flow pattern remains qualitatively similar to that 
of a medium made of rigid monodisperse particles. Sev¬ 
eral features observed in our model system have also 
been reported in other granular media composed of 


dry sand [34] and elastic disks |4|. The physical prop¬ 
erties of the individual particles used in these studies 
(size, shape. Young’s modulus and friction coefficient) 
are very different from the steel balls used in our ex¬ 
periments. These observations, taken in total, suggest 
universality of the observed flow patterns in granular 
media. 

The formation of a stagnation zone in flow past an 
obstacle has been reported for smaller grains [3411381 
[37] . In these cases they appear to occur only because 
of the existence of long-lasting interparticle contacts 
in a narrow region surrounding the obstacle. In our 
configuration, interparticle contacts are largely main¬ 
tained throughout the flow, because the particles are 
not in free fall. The occurence of the stagnation zone 
and its subsequent evolution (c/. Fig. |4(a)[ ) are deter¬ 
mined purely by interparticle rigid-body interactions 
and Coulomb friction at the contacts. The coefficient 
of friction determines the final extent of this stagnation 
zone, as confirmed by the close agreement between ex¬ 
perimental results and NSCD simulations (Fig. |4(b)| ). 
The occurence of a stagnation zone has implications 
for powder processing of metals and ceramics. In these 
applications, homogeneous densities are desired. In con¬ 
trast, in pharmaceutical powder compaction for making 
tablets, a gradation in density is sought. The stagna¬ 
tion zone inhibits consolidated particle mixing during 
compression and is undesirable both for uniform and 
graded densities. Thus, methods to reduce its size must 
be employed. Our experiments suggest that this may be 
achieved by reducing the coefficient of friction between 
individual grains (see Fig. |4(b)| ). This will reduce the 
base angle 0, and hence the extent of the stagnation 
zone, enabling free mixing. 

Classical Mohr-Coulomb plasticity solutions for quasi¬ 
static deformation also predict a stagnation zone [5] be¬ 
neath the punch, just prior to the onset of slip. The base 
angle of this zone is a function of the internal friction 
angle ip, and is 90° for materials with '0 = 0. While the 
present experiments and simulations have dealt with 
particle-level properties (such as the friction coefficient 
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/i), there is as yet no certain way to relate these to 
continuum-scale material constants like -0 m- 

Vortex-like circulation centers developed during the 
course of the flow (Fig. have also been reported in 
gravity-driven granular flow past obstacles [39] and as 
‘eddies’ in two-dimensional sheared granular materials 
mnn]. This is again a consequence of simple rigid body 
dynamics and friction, as seen from the NSCD simula¬ 
tions (c/. Figs. |5(a) and 5(b)[ ). Interestingly, as shown 
in Fig. vorticity in the material decays quickly with 
distance from the punch, reminiscent of a boundary 
layer in fluid flow. Such boundary layers have been re¬ 
ported in simulations of 2D simple shear m suggesting 
a correlation with large local dissipation due to inter¬ 
particle friction. These vortices may be reduced by the 
introduction of rounded edges on the punch, an idea 
that we hope to explore in the future. This could lead 
to the development of ‘laminar’ granular flows with re¬ 
duced energy dissipation, and better control of powder 
densities for processing applications. 

Shear band formation is common to many granu¬ 
lar materials containing different types of grains [35] . 
In our system, since the particles themselves were per¬ 
fectly rigid, plastic strain in the material resulted from 
large relative motion of neighboring particles. Regions 
of high velocity gradients also correspond to regions of 
high shear strain rate — see Fig.[^ As a consequence, if 
the kinematic configuration results in large velocity gra¬ 


dients (such as those marked by the arrows in Fig. 7(a)), 
the shear strain will tend to get localized in bands. In¬ 
side the band themselves, an increase in volume (dila- 
tancy) can be expected and this has been reproduced in 
the NSCD simulations (Fig.|^. However, in the present 
experiments, these bands developed as a consequence of 
competition between compaction due to the punch and 
material dilatancy |44| — while compaction of material 
immediately below the punch increases the local den¬ 
sity, the particles themselves prefer to move apart and 
create free space. 


Shear banding is a common cause of failure for Bulk 
Metallic Glasses (BMGs). The structure of many BMGs 
has often been approximated by hard sphere media na 


EHni, and the occurence of shear bands has been at¬ 
tributed to excess free volume concentration or the col¬ 
lective formation of STZs [8|. This idea is consistent 
with the present observations, since existence of local 
pockets of free volume can result in particle jumps, 
thereby inducing a large velocity gradient. 

The visual similarity between the vertical column 
underneath the punch (Fig. [IT] ) and a buckled elastic 
column shows how additional strain is actually accomo¬ 
dated. It is interesting to note that force chains in the 
medium also show a tendency to buckle [45] and most 
likely have a role to play in this process. One impor¬ 
tant distinction, however, is the reversibility in elastic 
column buckling — the deformation in Fig. [^is com¬ 
pletely irreversible. While further flow in the medium 
occurs after the bulge formation, this change cannot be 
reversed, even in the absence of gravity. 

Granular flows, both in the present (mesoscopic) 
case and fine-grained media (e.g., sand) [34], show re¬ 
markable similarity with punch indentation of rigid- 
plastic metals. In metal indentation, a central stagna¬ 
tion zone is developed right underneath the indenter, 
along with shear bands bounding it - exactly as in the 
present experiments. These have been predicted in met¬ 
als by slip line field analyses - the Prandtl field m - 
and confirmed in experiments [20]. Likewise, this flow 
is reproduced in exact-scale NSCD simulations of the 
granular medium. As in metals, the shear bands repre¬ 
sent regions of high local shear strain and significant ve¬ 
locity jumps (Fig.|^. The NSCD simulations also show 
that the formation of new surfaces when the punch pen¬ 
etrates the granular medium (Fig.[^ is very similar to 
that in metal indentation [36] . 

Given this strong similarity between indentation flows 
in metals and granular media, it is likely that the vortex- 
like structures observed in the present study also occur 
with metal indentation. While this is yet to be verified, 
perhaps, due to lack of adequate observations along in- 
denter walls in metals, it is nevertheless promising that 
such circulation zones have been captured in sliding 
metal surfaces [46] . 
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5 Conclusions 

The in situ imaging and 2D image analysis of punch 
indentation of a mesoscopic, dense granular medium 
has allowed flow features to be studied at high res¬ 
olution. Detailed information about particle positions 
and velocities has been obtained. The use of an anal¬ 
ysis method for the disordered grain network has en¬ 
abled the calculation of strain fields, velocity gradi¬ 
ents and rotations. Based on this data, the formation 
of collective patterns in the flow has been captured. 
Despite its coarse discrete nature, the model granu¬ 
lar material shows a number of flow features that bear 
striking resemblance to those seen in metals and fine¬ 
grained granular materials. A triangle-shaped stagna¬ 
tion zone is established underneath the punch, as in in¬ 
dentation of metals and sand. Vortices and shear bands 
are observed, and their characteristics have been quan¬ 
tified using experimentally measured deformation rates. 
These features are quantitatively reproduced using full- 
scale, 1:1 Non-Smooth Contact Dynamics (NSCD) sim¬ 
ulation, showing also that they originately purely from 
local rigid body dynamics and Coulomb friction. The 
coupled simulation has provided clues to the kinematic 
origin of shear bands, by direct visualization of the di- 
latancy in the plane of the bands. New flow features, 
such as free surface evolution and buckling of a vertical 
column of particles underneath the punch, not easily ac¬ 
cessible in experiments, are revealed by the NSCD. The 
observations in total suggest that, despite the large size 
of the particles, their collective dynamics very much re¬ 
semble a continuum. This is perhaps why continuum 
theories have been quite successful in analyzing quasi¬ 
static deformation problems pertaining to granular me¬ 
dia. 
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